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Abstract 

We investigate here the abiUty of a Green-Naghdi model to reproduce 
strongly nonlinear and dispersive wave propagation. We test in particular 
the behavior of the new hybrid finite-volume and finite-difference split- 
ting approach recently developed by the authors and collaborators on the 
challenging benchmark of waves propagating over a submerged bar. Such 
a configuration requires a model with very good dispersive properties, 
because of the high-order harmonics generated by topography-induced 
nonlinear interactions. We thus depart from the aforementioned work 
and choose to use a new Green-Naghdi system with improved frequency 
dispersion characteristics. The absence of dry areas also allows us to im- 
prove the treatment of the hyperbolic part of the equations. This leads 
to very satisfying results for the demanding benchmarks under consider- 
ation. 

1 Introduction 

The propagation of surface waves in an incompressible, homogeneous, invis- 
cid fluid is governed by the free surface Euler equations. In shallow water 
(i.e. when the typical wavelength of the waves is much larger than the typical 
depth), a good alternative is furnished by the Green-Naghdi equations (also 
called Serre or fully nonlinear Boussinesq equations); these equations were first 
derived in [5ni[IH[ini (see also [HI HI]), and we refer to [TS] for a recent review 
and a derivation under the formulation used here. It is known [1] that they 
approximate the full Euler equations with a good accuracy up to the breaking 
point. 

Let ho denote the reference depth, C(t, x) the elevation of the free surface 
with respect to its rest state, and b{x) the variation of the bottom topography 
with respect to the constant level z = —Hq. The Green-Naghdi equations couple 



the evolution of the water depth h{t, x) = ({t, x) + ho — b{x) to the evolution 
of u(t, x), the vertically averaged horizontal component of the velocity. For Id 
surface waves, these equations can be written 

(dth + d^{hu) =0, 

\dtihu) + d^ihu^) + ^ghd^C + [l + ahTj^r^i^ghd^C + /iQ(u)] = 0, 

(1-1) 

where a > 1 and the operators T and Q arc explicitly given by 

Tw = ~Y^lw - hd^hd^w + [d^Cdxh + ^dlb)w, (1.2) 

and 

Q{u) = 2hd^{h+^){d^uf+^h'^d^udlu+hdlhud^u+{d^C,dlh+^dlh)u^. (1.3) 

Remark 1.1. The usual Green-Naghdi equation corresponds to a — 1 in 
Using the classical observation that 

dtu = —dxC ~ udxU + higher order terms 

= adtu — (1 — a){dxC + udxu) + higher order terms, 

one deduces il.l}) . with the same accuracy (see for details). 

Looking at the linearization of (jl.ip around the steady state h = ho, u = 
over a flat bottom 6 = 0, one derives the dispersion relation associated to ()l.ip . 
It is found by looking for plane wave solutions of the form (/i, /m)e''-'^^~"'^ to 
the linearized equations, and is given by, 

2/,^ u ,2l + ("-l)(^^o)V3 
^Uk)=ghoe ,^^(,,^),^3 ■ (1.4) 

The interest of the parameter a introduced in Remark 11.11 is that it can be 
chosen in such a way that (jl.4p adequately matches the exact dispersion relation 
tt'^(fc) = gktanh{khQ) of the full Euler equation, even for non-small values of 
kho- Taking a = 1.159 as in [5] yields a good agreement up to kho = 4 on the 
linear phase velocity and up to kho = 2.5 on the linear group velocity. 

The goal of [S] is to develop a numerical code having good dispersive proper- 
ties, and able to handle successfully wave breaking and vanishing depth (shore- 
line) . The idea is to use a splitting scheme that decomposes the dispersive and 
nonlinear parts of the equations. More precisely, using a second order splitting 
scheme, the approximation J/""*"^ = (C""'"-'^, u""*"^) at time {n + l)St is found in 
terms of the approximation J7" at time n6t by solving 

C/"+i = 5i(<5t/2)52(<50^i('5t/2)f/", 

where S'i(-) is the solution operator associated to the nonlinear shallow water 
(or Saint- Venant) equations and S2{-) the solution operator associated to the 
dispersive effects. It is shown in [5] how to adapt this method to deal with 
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wave breaking and vanishing depth. Vahdations relying on analytical solutions 
and comparisons with experimental data give very good results. For most 
of the applications, the dispersive properties of this model (with a ~ 1.159 
in p.ip ) arc good enough. In order to consider challenging configurations 
that include high harmonics (typically up to khg = 4), we first derive a new 
family of Green-Naghdi models with improved frequency dispersion following 
the steps of Nwogu and adapt the approach of [5] to this new model. Since 
this is not relevant for this test case, we do not mention here the treatment 
of wave breaking nor of the shoreline (see [5] for details); on the other hand 
we propose numerical improvements on the computation of the hyperbolic part 
made possible by the absence of dry areas. 

The paper is organized as follows. In Section [2l we derive a new family 
of Green-Naghdi equations that have the same precision as the original one as 
kho — >■ 0, but with better dispersive properties. We then present in Section [3] 
the main features of the numerical scheme. Finally, we present in Section 3] 
some numerical validations of our derived models. 

2 A three- parameters family of GN models 

The main goal of this section is to derive a new family of Green-Naghdi equa- 
tions that are equivalent to the original set, up to higher order terms. As a 
mean to clearly identify these higher order terms in our fully nonlinear and 
weakly dispersive framework, we rather work on the nondimensionalized ver- 
sion of namely, 

(dth + d,ihu) = 0, 

I il + ^iahTj^)[dtihu) + ed^ihu") + ^hd^C] +ihd,C + e^ihQ{u) = 0, 

(2.1) 

where e = a/ho = 0{1) is the nonlinearity parameter and /i = h^/X^ ^ 1 is 
the shallowness parameter, a being the typical wave and bottom deformation 
amplitude, A the typical wavelength and ho the reference depth. We still denote 
by h the dimensionless water depth, h = 1 + e{( — b). 

In [18j . Nwogu showed that it was possible to improve the dispersive prop- 
erties of Boussinesq models by working with the velocity at a certain depth as 
dependent variable. This approach was generalized in [21] to the fully nonlinear 
case (Green-Naghdi). When the bottom is not flat, it turns out that in the fully 
nonlinear case, the standard Green-Naghdi equations written with the average 
velocity u do not belong to this new class of fully nonlinear models. This is the 
reason why we use a slightly different approach here, with the introduction of 
a new dependent variable ug that is not the velocity at a certain depth. The 
interest is that the computations arc somehow simpler and, more importantly, 
that the average velocity u appears as a particular case (9 = 0) — our choice of 
dependent variables also differs from the one used in jl6| and whose purpose is 
not to improve the dispersive properties but to work with potential variables. 
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For all > 0, wc thus define 

ug ^ {1 + fiery^u. (2.2) 

The first equation of (|2.ip can thus be written 

dth + d^{hug) + fxed^XhTue) = 0. (2.3) 

Note that this equation, as well as the first equation of p.ip is exact (and not 
up to higher order terms). For the second equation, the transformation is more 
complex. Let us first focus on the component (1 + iiahTj^)dt{hu). Using (|2.2p 
we can write 

dtihu) = {1 + ^i9hT^)dtihue) + fie[dt,hT^]hue. (2.4) 

Using (jl.2|) . we deduce after some computations that 

12 1 

[dt,hTj]w = --hdthdlw - -dx{hdth)dxw 

+ {^dl{hdth) + dx{dth)dxb + ]^{dth)dlb)w- 
with the help of p.3|) . this shows that 1(2. A\ can be put under the form 

dt{hu) = (1 + fiehT^)dt{hue) + en0Qi{hue) + O(m^), (2.5) 

with 

Qi(m) = ^hdx{ml) + ^{dxh)ml~^mdl{hmx)~{mxxdxb+^mxd'^b)m. (2.6) 

Let us now focus on the component e(l + fiahTj^)dx{hu^) . Using (|2.2p . one 
has 

edx{hv?) = edx{hul) + eisOdxi^huoTue) + 0{^i^), 

= e(l + ^i6hT^)dx{hul) + EfieQriug) + 0{fi^), 

with Qr = dx{h{2ugTu0 — TuO'^)) + [dx, hTj^]{hul). After some computations, 
we obtain Q, (we) = hQ{ug) + Q^iug) with 

Q2{ug) = -^dx{h^)dl{hul)-(^dl{h')^edxQdxb)dx{hul)~h''dxb{dxUef 
+ (^dl{h^) + dx{2edxCdxb+^dlb) - dxbdlbyul (2.7) 
We have thus proved that 
edxihu") = e[l + ^i9hr^)dx{h4) + efie{hQ{ug) + Q2{ug)) + 0{^l^). (2.8) 
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Pluging (|2.5p . (|2.8p into the second equation of (|2.ip and recalling (|2.3p . we thus 
get that (|2.ip is equivalent, up to terms in the second equation, to 

dth + edx{hug) + e^9dx{hTug) — 0, 

(1 + fi{a + 9)hTj;)[dt{hue) + edx{hul) + ^hd^Q] (2.9) 
-^T^hd.C + etJi{{l + 9)hQ{ug) + BiQiihue) + Q2iue))) = 0. 

At this step, we point out that the introduction of the parameter 9 reveals to be 
useless as far as the minimization of the dispersion errors is concerned. Indeed, 
keeping in mind that a > 1 and 9 > 0, computing the dispersion relation 
associated to these new equations and minimizing the dispersion errors (on 
the phase and group velocities) over the range kh G [0,4] yields the values 
a = 1.159 and 9 = 0. 

A workaround for this issue is to introduce a third parameter 7 > by applying 
the operator {1 + fi-yhTj-J to the first equation of (|2.9p . Doing so, the first non- 
dimensional equation becomes 

dth + ed^ihug) + efi9{l + n-ihT\y^dx{hTug) = 0{fi^). 

h 

The drawback of this manipulation is that the exactness of the equation is lost; 
nevertheless, the numerical computations show that it is a fair price to pay for 
the considerable improvement of the dispersive properties thus achieved. 
On the whole, this final manipulation yields a new three-parameter family (Q) 
of Green-Naghdi models Qa,e,'y that can be written in dimensional form as 



dth + d.,{hug) + 9{1 + ihT\)-^dx{hTug) = 0, 

h 

dtihug) + dAhul) + ^^ghdxC (2-10) 
a + 9 
.l._i/l " 



+ (1 + {a + 9)hT-)-^ (^niQhdxC + hQ{ug)) = 0, 
h \a + 9 J 

with a > 1, 6* > and 7 > and 

Q(ue) = (1 + Q)Q,{ug) + ^-{Qx(]-iug) + Qiiug)). (2.11) 

Remark 2.1. We insist on the fact that all these models are equivalent at 
order 0{iJ?) to the original Green-Naghdi equations (OUP, which correspond to 
the particular case a = 1, 9 = ^ = Q. 

Now, one easily computes the dispersion relation associated to (|2.10p . 

Starting from this expression, one can deduce the expression - not reported 
here - of the associated linear phase and group velocities, denoted C^jy and 
Cq^, along with the linear shoaling coefficient 7^. Using these expressions, one 
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computes for each triplet (a, 9, 7) the corresponding dispersive errors between 
these model properties and the theoretical ones coming from Stokes theory, 
namely C^, Cg and 70 : see [SIE] for further details. The minimization of the 
shoaling error is quite problematic as a better shoaling involves a deterioration 
of the phase and group velocities, and vice-versa. For this reason, we provide 
here two different sets of optimized triplet (a, 0,7), one corresponding to op- 
timal phase and group velocities errors only - which can be seen as optimal 
values for flat bottoms - and one corresponding to optimal phase and group 
velocities and shoaling errors. Optimizing over the dispersive range fc/io £ [0, 4] 
yields : 

dopt = 1.028, 6lopt = 0.188, 7opt = 0.112, for flat bottoms, (2.12) 

and 

a„^pt = l, 6l„^pt = 0.207, 7„''pt = 0.071, for uneven bottoms. (2.13) 

As shown in Figure [U the dispersive range of validity of the model cor- 
responding to the first set is considerably extended (up to /c/iq = 8 for the 
phase velocity and up to fc/iQ = 5 for the group velocity) but with deteriorated 
shoaling properties (satisfactory up to fcft-o = 1.5 only). On the other hand, 
the range of vahdity of {a opt, 6 opt :1 opt) model is similar to the one of the one- 
parameter model if we consider the phase and group velocities, but exhibits a 
very good shoaling up to kho = 4. 



3 Numerical scheme 

We present here the main lines of our numerical scheme. The splitting scheme 
is explained in ^13. 11 we then show in Wi.2\ and i j3.3l how we treat respectively 
the hyperbolic and dispersive parts of the equations. 



3.1 The splitting scheme 

We decompose the solution operator S{-) associated to (|2.10p at each time step 
by the second order splitting scheme 

S{6t) = Si{St/2)S2{St)Si{St/2), (3.1) 

where 5*1 and S2 arc respectively associated to the hyperbolic and dispersive 
parts of the Grcen-Naghdi equations (|l.ip . More precisely: 

• Si (t) is the solution operator associated to the nonlinear shallow water equa- 
tions 

dth + dxihu) ^0, , , 

dt {hu) + {hu^ + gh^/2) = -ghd^^b. ^ ' 

• S2{t) is the solution operator associated to the remaining (dispersive) part 
of the equations, 

dth + eil + -fhTj^)-^d^ihTug)^0, 

dt{hu) - ^ff/^a.C + (l + (a + 0)/iri)-i [^ghd^C + hQiue)] = 0. 

(3.3) 

where the operators T and Q are as in (|1.2|) and ()2.11|) . 
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Figure 1: Errors on linear phase (top) and group velocities (bottom), and 
shoaling coefficient. The reference from Stokes theory is in black plain line, 
the standard GN model (a = 1) is in dots, the optimized one-parameter model 
(a = 1.159) in dash-dots, the three-parameter model optimized for flat bottom 
p.l2p in grey plain line, and the three-parameter model optimized for uneven 
bottom (|2.13|) in dashes. 



As far as the time discretization is concerned, we use a fourth-order explicit 
Runge-Kutta time scheme, both for the hyperbolic and the dispersive part. 

3.2 Spatial discretization of the hyperbolic part S'i( ) 

We use a high order finite volume approach in conservative variables w = 
* {h,hu), relying on Riemann problems for hyperbolic conservative laws |11| . 
This allows accurate computation of propagating bores, with reduced spurious 
effects of numerical dissipation and dispersion (this property is used in [5] to 
handle wave breaking). Since we aim at computing the complex interactions 
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between propagating waves and topography (including the preservation of mo- 
tionless steady states), we embed this approach into a well-balanced scheme. 

Based on discrete finite- volume cell averaging at time t" = n6t, we use 
5*''-order accuracy WENO reconstructions, following [T3], together with the 
weight splitting method [13]. Considering a cell Ci, this approach provides, 
for all i", interpolated quantities Wi_/, w^^c and w^ ,., respectively at the left 
boundary, center and right boundary of the cell. To get a stable and well- 
balanced scheme, the following reconstructions are introduced, following [5]: 

b* = max(fei,r,&j+i,z), 

h*r = max(0, h,^r + bi^r ~b*), 

K+i.i = max(0, /li+i,/ + - b*). 

These left and right values for h* are used to provide auxiliary values w* ^ and 

which are injected into a Riemann solver. The corresponding semi-discrete 
finite-volume scheme for p.2[) reads: 

^W,(t) + ^(|'I'''(wi,r,Wi+i,;,6i,r,&j+l,;) - {'Vfi-l.r ,^i,ubi-l.r M,l)^ = Sc,i 

where Vl/'' and are numerical flux functions based both on a conservative 
flux consistent with the homogeneous SWE and the correction to the interface 
fluxes due to the hydrostatic reconstructions. Scs is a centered discretization 
of the source term needed to achieve consistency and well-balancing properties. 
The reader is referred to [5] for any details concerning the hydrostatic recon- 
struction and to [3] for details concerning the conservative flux, issued from 
a relaxation approach, used in our splitting method. To achieve an overall 
4*'' order accuracy, following [17], we define 5'c,i = *(0,S'^"), with use of the 
following quadrature rule |17j : 

Scl = |(4((/j«J-f/j,.c)(fe,J-6.,c) + (/l^,c + /ti,i)(^«,c-&j,0)-(/li,i+/i»,r)(6»,/-&«,))- 

This discretization allows both high-order accuracy and well-balancing prop- 
erties. If the benefits are obvious concerning accuracy and convergence rates, 
the main drawback is the lack of robustness of this approach. In particular, we 
have to avoid any situation in which dry areas can occur, as in [5]. 

3.3 Spatial discretization of the dispersive part S2{-) 

As specified in the system (|3.3p is solved at each time step using classical 
fourth-order accuracy finite-differences. For the sake of clarity, the reader is 
referred to [5] for details concerning discretization of these dispersive terms, 
and also for the important issue of boundary conditions. 

The choice of a finite difference method for solving 5*2 also entails to switch 
between the cell-averaged values and the nodal values of each unknown, in a 
suitable way that preserves the global spatial order of the scheme. 
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4 Numerical validations 



4.1 Propagation of nonlinear cnoidal waves 



It has recently be shown in [6] that system (|2.10p with a = l,6' = 7 = admits 
the following family of periodic cnoidal solutions: 



h{x,t) = flo + aidn^(K(a; — ct), fc), 
ho 



(4.1a) 
(4.1b) 



/3ai 



21/^0(00 + ai)(ao + (1 - k'^)ai) 



y/gaoiao + ai)(ao + (1 - fc^)ai) 
ho 



where k G [0, 1], ag > 0, ai > are real parameters and dn(-,fc) is a Jacobi 
elliptic function with elliptic modulus k. The parameters of this solution can 
be related to physical variables in order to define (|4.ip in terms of wave height 
H, wave period T, and mean water depth ho. This can be achieved by solving 
the equations: 



ai 



H 

fc2' 



ao^ ho - ai 



E{k) 

K{ky 



4[aoK{k) + aiE{k)Y 



(4.2a) 



where cij = 2Tr/T is the angular frequency, K{k) and E{k) are the complete 
elliptic integrals of the first and second kinds respectively. 
In this test case, we study the propagation of strongly non-linear cnoidal 




Figure 2: Propagation of a periodic cnoidal wave over a periodic domain {H = 
0.6 TO, /iQ = Ito, T = 4s). Exact solution is plotcd in blue line, numrical 
results at i = 60 s in red line. 
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waves defined with H = 0.6 m, /iq = Im and T = As. The computational 
domain length is equal to 2 wave-lengths and we use periodic boundary con- 
ditions. In Figure ^ numerical and theoretical solutions are compared after 
15 periods of propagation, at t ^ 60s. Computation is performed with Sx = 
and a Courant number equal to 1. At t = 60s, we obtain a relative amplitude 
error of 1.3 10~^ % and the relative celerity error is estimated to be less than 
1.010~^%. Very accurate results are thus obtained, as an assessment of the 
ability of our numerical method to compute in a stable way the propagation of 
non-linear cnoidal waves. 

4.2 Propagation of highly dispersive waves 

x(m) 

2 4 6 8 10 12 14 16 18 20 22 24 




Figure 3: Experimental set-up and locations of the wave gauges 

In this section, we investigate the ability of the previous models to describe 
the propagation and the interaction of highly dispersive waves. To this end, 
we examine the propagation of regular waves over a submerged bar, using 
the set-up introduced by Beji & Battjes (1993, [3]), and first used as a test 
by Dingemans (1994, ^Q\), in which the author evaluated the performance 
of various Boussinesq-type formulations by comparing computed free-surface 
time-series with experimental measurements at several gauges (see Figure [3]). 

Wc consider here Case A of [TU] where a small-amplitude long wave is 
generated at the left boundary : the free-surface elevation a is 0.01m, the 
time period T is 2.02s and the initial depth /ig is 0.4m, which corresponds 
to a dispersion parameter kh^ « 0.7. When the incident wave encounters 
the upward part of the bar, it shoals and steepens, which generates higher- 
harmonics as the nonlinearity increases. These higher-harmonics are then freely 
released on the downward slope, and become deep-water waves behind the bar. 
As discussed in Woo & Liu (2001, [22]), significant wave energy is present at 
kho « 4 in the region behind the bar. For this reason, models based on a 
weakly dispersive assumption - such as Boussinesq-type ones or Green-Naghdi 
equations - are generally not able to correctly reproduce the measured profiles, 
since their linear dispersion properties usually become inaccurate beyond kho w 
3. As it involves highly nonlinear interactions and requires extended dispersive 
properties, this discriminating test has become a widely used benchmark over 
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the past decade. 

Numerical simulations for this case are shown in Figure |4l where computed 
free-surface time-series are compared to experimental measurements at the last 
four gauges, located on the downward part (#8) of the bar and behind it (#9, 
#10 and #11). Results at gauges #1 to #7 are not reported here as they do not 
exhibit any significant difference between computations and experimental data. 
We choose here to compare the numerical results of two Green-Naghdi models: 
the optimized one-parameter model (|l.ip with a = 1.159 and the optimized 
three-parameter model (|2.10p with a = 1, 6* = 0.207 and 7 = 0.071. We point 
out that the optimized three-parameter model for flat bottoms (|2.12p yields 
results that are very similar to the optimized one-parameter model, and thus 
not reported here. For these simulations, we used the numerical parameters 
6x = 0.03m and St = 0.01s, which corresponds to a Courant number equal to 
0.6. The incident wave is generated at the left boundary using an absorbing- 
generating method presented in [S], and a sponge-layer is used at the right 
boundary to prevent perturbations by reflected waves. Computed time-series 
are shown over the time window [35s, 39s] where the waves kinematics are well 
established. 

Both models show a similar and very satisfactory agreement with experi- 
mental measurements on the gauge #8 located on the downward part of the 
bar, which is an expected result as the higher harmonics are not released yet. 
The first discrepancies between the two models appear on gauge #9 located 
after the bar, where the three-parameter model provides a slightly better agree- 
ment with the measured profile than the one-parameter model. This remark 
persists on gauge #9, even if both results still perform well. However, a very 
significant difference appear on the computed profiles at the very last gauge : 
the one-parameter model does not reproduce the correct profile (the amplitudes 
are either under-estimated or over-estimated) while the three-parameter profile 
shows an overall very good agreement with the results. This last gauge is by far 
the most discriminating one as the higher-harmonics are completely released at 
this location and can be seen as deep-water (and thus highly dispersive) waves. 
Keeping this in mind, the results obtained with our three-parameter model ap- 
pears to be an excellent performance and clearly shows that the optimization 
presented in (|2.13[) was very helpful in reproducing the correct profile. More- 
over, seeing that the three-parameter model optimized for flat bottoms (|2.12p 
yields less accurate results, even if this model exhibits significantly better dis- 
persive properties on the linear and phase velocities, we can reasonably assume 
that the linear shoaling optimization is an unavoidable procedure for this very 
challenging test case. 
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One-parameter optimized model Three-parameter optimized model 




35 36 37 38 39 35 36 37 38 39 



t (s) t (s) 

Figure 4: Comparison of computed free-surface time-series with experimen- 
tal measurements for the optimized one-parameter Green-Naghdi model (left) 
and the three-parameter Green-Naghdi model optimized for uneven bottoms 
(right). 
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